Calhoun:  The  NPS  Institutional  Archive 
DSpace  Repository 


Reports  and  Technical  Reports 


All  Technical  Reports  Collection 


1994-06 

Efficient  parabolic  equation  solution  of 
radiowave  propagation  in  an  inhomogeneous 
atmosphere  and  over  irregular  terrain  :  formulation 

Janaswamy,  Ramakrishna 

Monterey,  California.  Naval  Postgraduate  School 


http://hdl.handle.net/10945/28720 


Downloaded  from  NPS  Archive:  Calhoun 


DUDLEY 

KNOX 

LI&RARY 


CalliDun  is  a  proiect  of  the  Dudley  Knox  Librarv  -at  NPSj  furthering  the  precepts  ?nd 
gioali  of  open  government  ai^d  govemmenl  tFansparencv.  All  inlbrmatfoni  contained 
herein  has  been  approved  for  release  by  the  NPS  Public  Affairs  Officer. 

Dudtey  KnoK  Ubrary  /  Naval  Postgraduale  School 
411  Dyer  Road  /  1  Univefsity  Circle 
MontereYf  California  USA  93943 


http://w  w  w,n  ps.edu/litiTary 


NPS-EC-94-004 


NAVAL  POSTGRAOOATE  SCHOOL 

Monterey,  California 


Efficient  Parabolic  Equation  Solution 
of  Radiowave  Propagation  in  an 
Inhomogeneous  Atmosphere  and  Over 
Irregular  Terrain:  Formulation 

by 

Ramakrishna  J  anaswamy 
June  1994 


Approved  for  public  reslease;  distribution  unbmited. 
Prepared  for:  Naval  Postgraduate  School,  Monterey,  CA 

FedDocs 
D  208.14/2 
NPS-EC-94-004 


Naval  Postgraduate  School 

Monterey,  California  93943-5000 

Rear  Admiral  T.  A.  Mercer 
Superintendent 

This  report  was  prepared  for  the  Naval  Postgraduate  School. 
Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  report  was  prepared  by: 


H.  Shull 
Provost 


Department  of  Electrical  and 
Computer  Engineering 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704^0168 


bvroen  tK»%  oi  fntor>n»Uon  lO  *^er»Q*  <  ^Owf  Ot'  ft\POn\f,  tht  t«me  lo»  mtiruritoni.  te*r<hin^  •vritm^  d«t«  »Ow»<ei. 

a#the»*»'0  th9  4*l*  »nd  (0*rtpletrn^  0f*d  tf>t  <0**tnr0*f  Of  *«lo»m*!»on  CO'V'm^nit  r*o«rd«'^  iKn  bvrd^n  Of  •«y  OlK^r  of  thM 

toi»ect*o«  ♦nlOfm*f»Of'.  •fHtudfog  iu99«-itk>fn  lof  fedo<if>9  IKn  bwfde"  to  w*tK*f»9ton  Mt»*dou«r!eri  Servfce^.  D'ffnof*te  »of  »nfofm*tiOf»  ODef«t«on%  •f»d  «*DOft».  WiS  Jefferson 
1 W.  AfliKqlon.  vn  2}m^y0}.  »»^X©  fKi*  0»«<f  O*  JKd  ►wrwort  Wfdoaioo  l»roiecl  (OrTJa  O  111).  >A/»$Kif*9tor>.  OC  /DSOJ 


1.  AGENCY  U$E  ONLY  (Lesve  bUnk) 


2.  REPORT  DATE 

June  1994 


3,  REPORT  TYPE  AND  DATES  COVERED 

Technical  Report  1  Jan  94-1  Jun  94 


4.  TITLE  AND  SUBTITLE 

Efficient  Parabolic  Equation  Solution  of  Radiowave  Propagation  in 
an  Inhomogeneous  Atmosphere  and  Over  Irregular  Terrain:  Formulation 


6.  AUTHOR(S) 

Ramakrishna  Janas wamy 


5.  FUNDING  NUMBERS 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS{ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

NPS-EC-94-004 


9.  SPONSORING/MONITORING  AGENCY  NAME{S)  AND  AOORESS{ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943 


10.  SPONSORING/MONITORING 
AGENCY  report  NUMBER 


11.  SUPPLEMENTARY  NOTES 

The  views  expressed  in  this  report  are  those  of  the  authors  and  do  not  reflect  the 
official  policy  or  position  of  the  Department  of  Defense  or  the  United  States  Government. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  public  release; 
distribution  is  unlimited. 


12b.  DISTRIBUTION  CODE 


13.  abstract  (M3Mtmum  200  word^ 

Formulation  is  given  for  efficient  parabolic  equation  solution  of  radiowave  propagation  in  inhomogeneous 
atmosphere  and  over  irregular  terrmn.  Both  standard  and  wide  angle  parabolic  equation  derivations  are 
presented.  Impedance  boundary  conditions  are  used  to  characterize  the  ground.  A  tropospheric  boundary 
condition  based  on  the  exact  solution  of  Schrodinger  equation  in  a  quarter  plane  is  derived.  To  permit 
efficient  modeling  of  the  irregular  boundary,  the  parabolic  equation  together  with  the  boundary  conditions 
is  transformed  into  a  numerically  generated  curvilinear  coordinate  system.  Finally,  formulation  is  presented 
for  a  finite  difference  solution  using  Crank-Nicolson  implicit  scheme. 


If.  SUW*CT  TEfWS  

Kadiowave  Propagation,  Parabolic  Equation,  Finite  Differences 

15.  NUMBER  OF  PAGES 

37 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

u^cYa'^ified 

18.  SECURITY  CLASSIFICATION 

19.  SECURITY  CLASSIFICATION 

20.  LIMITATION  OF  ABSTRACT 

SAR 

NSN  T540.01-280-5500 


Standard  Fo»ni  298  (Rpv  2-89) 

by  ansi  Sid 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOS 
MONffiREV  OA  03043  6101 - - 


/ 


Abstract:  Formulation  is  given  for  efficient  parabolic  equation  solution  of  radiowave 
propagation  in  inhomogeneous  atmosphere  and  over  irregular  terrain.  Both  standard 
and  wide  angle  parabolic  equation  derivations  are  presented.  Impedance  boundary 
conditions  are  used  to  characterize  the  ground.  A  tropospheric  boundary  condition 
based  on  the  exact  solution  of  Schrodinger  equation  in  a  quarter  plane  is  derived.  To 
permit  efficient  modeling  of  the  irregular  boundary,  the  parabolic  equation  together 
with  the  boundary  conditions  is  transformed  into  a  numerically  generated  curvilinear 
coordinate  system.  Finally,  formulation  is  presented  for  a  finite  difference  solution 
using  Crank-Nicolson  implicit  scheme. 
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1.  Introduction 


It  is  well  known  that  multipath  fading  can  significantly  effect  the  link  reliability  in  a 
communications  system  or  target  detectability  in  the  case  of  radar.  The  path  between 
a  transmitted  and  a  receiver  is  often  obstructed  by  natrual  or  man-made  obstacles 
such  as  hills,  buildings,  atmospheric  layers,  trees,  rain,  fog,  etc.  Propagation  outage 
due  to  multipath  fading  depends  in  a  complicated  manner  on  propagation  climate, 
terrain  features,  path  length,  radio  frequency,  and  fade  margin.  In  the  case  of  atmo¬ 
spheric  multipath  fading,  interference  due  to  two  or  more  super- refracted  rays  arriving 
at  the  receiver  via  different  paths  can  lead  to  a  complete  loss  of  signal.  Other  patho¬ 
logical  phenomenon  such  as  obstruction  fading  (caused  by  sub-refractive  atmospheric 
effects)  and  ducting  (caused  by  extreme  super-refractive  effects)  are  also  possible. 
Such  phenomenon  are  more  common  in  warm  and  tropical  climates,  particularly  near 
shores,  where  elevated  inversions  are  formed  easily  due  to  the  large  temperature  and 
partial  pressure  differentials.  Reflection  multipath  fading,  which  is  due  to  interfer¬ 
ence  between  the  direct  and  the  ground  reflected  ray  depends  strongly  on  the  terrain 
geometry  and  ground  constants.  Moreover,  elevated  terrain  features  could  completely 
mask  a  receiver  from  a  transmitter  leading  to  severe  loss  of  signal  (in  some  cases,  it 
is  advantageous  to  site  antennas  behind  hills  to  provide  shielding  against  undesirable 
interference).  It  is  very  important  to  assess  the  effects  of  environment  on  the  link.  A 
computer  model  that  can  take  into  account  a  given  refractive  index  profile,  terrain 
elevation  data,  and  varying  ground  parameters  will  be  very  helpful  in  predicting  the 
link  performance. 

In  this  report  we  present  formulation  details  for  an  efficient  numerical  solution  of 
wave  propagation  in  an  inhomogeneous  atmosphere  and  over  irregular  terrain  using 
parabolic  equation. 

Unlike  all  other  previous  formulations  of  the  parabolic  equation,  we  will  use  a  modified 
Helmholz  equation  for  propagation  in  an  inhomogeneous  atmosphere  as  suggested  by 
Maxwell’s  equations  (all  previous  formulations  use  a  Helmholtz  equation  which  is 
only  true  for  fields  in  a  homogeneous  medium).  Because  the  parabolic  equation  is  a 
full-wave  method,  it  will  include  all  aspects  of  wave  propagation  such  as  reflection, 
refraction,  diffraction,  and  surface  wave  propagation.  In  this  respect  it  is  far  superior 
to  the  commonly  used  ray  method. 

Parabolic  equation  approximation  to  an  elliptic  partial  differential  equation,  which 
the  true  fields  satisfy,  has  proven  to  be  a  viable  approach  for  studying  propagation 
problems  in  underwater  acoustics.  The  method  is  just  gaining  popularity  with  the 
electromagnetic  community.  Although  the  parabolic  equation  regards  waves  as  es¬ 
sentially  traveling  one-way,  it  allows  a  rapid  solution  of  the  fields  by  way  of  marching 
along  the  range  starting  from  an  initial  range.  Another  advantage  of  the  PE  method 
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compared  to  the  ray  methods  is  that  it  is  valid  even  in  the  shadow  region  where  the 
simple  ray  methods  completely  break  down.  Furthermore,  it  appears  to  be  the  only 
practical  method  for  predicting  propagation  over  long  ranges  (greater  than  1  km)  over 
a  wideband  from  HF  (a  few  MHz)  through  SHF  (a  few  tens  of  GHz).  The  method  is, 
however,  not  without  limitations.  In  its  standard  form,  the  accuracy  of  the  method 
is  limited  to  waves  traveling  essentially  within  ±10°  from  horizontal.  Furthermore, 
treatment  of  the  boundary  conditions  on  the  uneven  terrain  is  difficult. 

The  method  we  propose  to  use  will  attempt  to  remove  both  of  these  defficiencies. 
Firstly,  we  use  a  Helmholtz-like  elliptic  equation  to  describe  the  fields  and  arrive 
at  a  wide-angle  parabolic  equation  subject  to  certain  approximations.  To  facilitate 
imposition  of  the  boundary  conditions  on  the  irregular  terrain,  the  equations  will  be 
transformed  to  a  body-fitted  curvilinear  coordinate  system.  The  PE  will  be  solved 
using  finite  differences  on  a  non-rectangular  mesh.  This  is  a  major  departure  from 
previous  approaches  which  will  not  only  make  the  method  more  efficient  but  also 
more  accurate. 

In  Section  2,  we  derive  the  exact  equations  satisfied  by  2D  fields  in  an  inhomogeneous 
atmosphere.  Impedance  boundary  conditions  are  used  to  characterize  the  ground.  In 
Section  3,  we  present  details  on  the  impedance  boundary  conditions.  Starting  from 
the  exact  equations  presented  in  Section  2  for  the  fields,  we  derive,  in  Section  4,  a 
parabolic  equation  (PE)  valid  for  narrow  angle  propagation.  This  case  is  termed  as 
the  standard  PE.  The  standard  PE  is  generally  valid  for  propagation  angles  that 
are  within  ±10°  from  horizontal.  To  accomodate  waves  traveling  at  larger  angles,  we 
present  the  derivation  of  a  wide  angle  PE  in  Section  5.  To  truncate  the  computational 
domain  we  derive  boundary  conditions  on  an  upper  boundary,  which  are  termed  as 
the  tropospheric  boundary  conditions.  The  derivation  is  based  on  the  solution  of 
Schrodinger  type  parabolic  equation  in  a  quarter  plane  x  >  0,  i/  >  0.  This  is  presented 
in  Section  6. 

For  an  efficient  numerical  implementation,  we  transform  the  differential  equation  and 
boundary  conditions  to  a  curvilinear  coordinate  system.  This  is  presented  in  Section 
7.  Details  on  the  numerical  generation  of  the  curvilinear  coordinate  system  are  given 
in  Section  8.  Finally,  in  Section  9,  we  present  steps  leading  to  a  finite  difference 
solution  of  the  equations  using  a  Crank- Nicholson  implicit  scheme. 
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2.  Solution  of  2D  Fields  in  an  Inhomogeneous  Medium 


Consider  an  electric  source  producing  fields  in  an  inhomogeneous  region  as  shown 
in  the  figure  below.  Let  us  assume  that  both  the  sources  and  the  medium  are  two- 
dimensional  in  nature  in  that  all  quantities  are  independent  of  the  ^-coordinate.  As 
in  the  case  of  a  homogeneous  medium  the  fields  can  be  decomposed  into  slTEz  case 
(vertical  polarization)  and  a  TMz  case  (horizontal  polarization).  It  is  assumed  that 
propagation  takes  place  in  the  xr/-plane. 


From  Maxwell’s  equations,  we  have 

V  X  E  =  —jujfiH 


V  X  //  =  j^tE  +  J 


(1) 

(2) 


In  the  case  of  vertical  polarization,  the  fields  could  be  written  in  terms  of  the  2- 
component  of  the  magnetic  field,  H  =  i//z,  {T Ez  fields).  Substituting  into  (2)  we 
have 

V  X  H  =  V  X  (zHz)  =  VZ/z  X  z  —  ju)cE  +  /=>  cE  = - — - — — 


In  a  source- free  environment,  we  have 


cE  = 


VHz  X  ‘z 


Substituting  into  (1)  we  get 


V  X 


E  = 


ju 


=  —jiijfiH 
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The  equation  satisfied  by  the  magnetic  field  is  then 


— V  . 

jio  \  e 


^  =  -juiiHz 


V  •  (^)  +  =  0 


Vertical  Polarization 


Once  the  magnetic  field  is  determined  the  electric  field  is  given  by 


jue 


Note  that  does  not  satisfy  the  Helmholtz  equation  unless  c  is  constant. 
For  horizontal  polarization  on  a  similar  analysis  with  E  =  zE^  shows  that 


Horizontal  Polarization 


(3) 

(4) 

(5) 


If  the  medium  is  non-magnetic,  fi  =  fio  and  E^  satisfies  the  Helmholtz  equation.  We 
will  characterize  the  ground  in  terms  of  impedance  boundary  conditions  [3]. 

We  may  combine  the  vertical  and  horizontal  cases  shown  in  (3)  and  (5)  into  an 
equation  of  the  form 

V  •  (aVt/))  +  /?V’  =  0  (6) 

where 


-  TE  or  Vertical  Pol. 
e 


a  =  < 

—  TM  or  Horizontal  Pol. 

(7) 

fi 

^  =  akln^{x,y)>0 

(8) 

'  Hz  TE  Pol. 

'll;  =  < 

^  Ez  TM  Pol. 

(9) 

The  quantity  n(x,  y)  is  a  position  dependent  refractive  index  of  the  medium.  The 
partial  differential  equation  (PDE)  given  in  (6)  is  elliptic,  for  we  have 


d  (  d'ip\  d  (  dil’\ 

^  j  ^  V  ^ j 


+  /0  V’  =  0 
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or 


a 


34)  da  dip  da 
+  1-  +  /?0  =  0 
dx  dx  dy  dy 


dx'^  dy^ 

The  discriminant  for  the  above  PDE  is  —1  <0  implying  elliptic  nature  of  the  equation. 
Equation  (6)  could  be  expressed  as 


,  I  da  \  da  B 

^  ^  +  -  ~^'^y  +  -^  =  0 

dx  a  dy 


a 


a 


For  a  non-magnetic,  loss-less  medium,  we  have 

1 


o  =  < 


Vertical  Pol. 


— ,  Horizontal  Pol. 


CqC 

1 


so  that  for  vertical  polarization 


1  da 
a  dx 


d  (\ 


^  dx\tr)  tr  dx  11?  dx 

2  dn  .... 

= - —  [n  ]s  the  refractive  index) 

n  dx 


Similarly, 

Letting 


1  da 
a  dy 


2  dn 
n  dy 


ai(x,y)  =  < 


a2{x,y)  =  < 


2  dn 
n  dx 

0 

2  ^ 

n  dy 

0 


Vertical  Pol. 
Horizontal  Pol. 
Vertical  Pol. 

Horizontal  Pol. 


equation  (10)  may  be  expressed  in  the  form 

V^i/’  d"  Oi0i  +  (i2ipy  +  k^n^ Ip  =  0 


(10) 


(11) 
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3.  Impedance  Boundary  Condition 


Impedance  boundary  condition  relates  the  tangential  components  of  electric  and 
magnetic  fields  at  the  interface  of  two  media.  If  i)  is  a  unit  normal  and  i  is  a  unit 
tangent  as  shown  in  the  figure  below,  the  boundary  conditions  are  given  by  [3] 


A 

s 


s 


s  =  X  cos  9  y  sin  6 
V  =  —X  sin  9  -\-  y  cos  9 

X  =  s  cos  9  —  0  sin  9 
y  =  s  cos  9  0  cos  9 

I  (12) 

where  =  Zsfyo  is  the  surface  impedance  normalized  to  the  free  space  value  7]o. 
The  equation  may  also  be  expressed  as 


(0  •  E)0  —  E  =  — r/oAjz)  x  H  =>  0  x  E  =  TjoAgO  x  (0  x  H) 

=  7/oAs  [(i>  ■  H)if  -  a 


(13) 


The  surface  impedance  is  determined  from  the  intrinsic  impedance  of  the  medium 
by  considering  plane  wave  reflections  from  the  interface.  The  complex  propagation 
constants,  71,  72,  and  the  intrinsic  impedances  t}i  and  72  in  terms  of  the  media 
constants  are  indicated  in  the  figure. 


7?  = 


ju!yrl^lo{crl  +jAcoeri) 
-kofiri  (en  -i— ) 

^of^r2^Tc2 


According  to  Snell’s  law,  we  have  71  sin  9,  =  72  sin  9t 
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The  plane  wave  reflection  coefficients  for  the  vertical  and  horizontal  polarizations,  Ry 
and  R}j  are  [7] 

7/2  sec  61  —  r/i  sec  6^ 


Rh 

Ry 


r/2  sec  6t  +  rji  sec  62 

7/2  cos  9t  —  7]i  cos  Ot 
rj2  cos  6t  —  Tji  cos  6^ 


Surface  Impedance  =  ri2  sec  6t 
=  t]2  cos  6t 


= 


7]2  sec  Ot  7]2 

m  vi 


1  -  |^^sin0, 


-1/2 


/  I^t2^tc\ 

.  /^rl  ^rcl  2  / 

1  COS  rpi 

J 

y  /^rl  ^tc2 

/^r2^rc2 

and 


A''  = 


1-1/2 


/  f^T2^TC\ 

1 - COS 

J 

y  ^rc2 

I^t2^tc2 

1/2 


For  the  special  case  of  pi\  =  H2  =  f^o-  <^i  =  0, 

1  - 


Af  = 


Ar  = 


^r2  —  J^t2 

For  normal  incidence  t/j,  =  90° 


1-1/2 


COS^  Xpi 


1  1/2 


cos^ 


(14) 

(15) 


Af  =  A^  =  J 

V  £r2  -jcrr2 

For  the  2-D  case,  impedance  boundary  conditions  for  vertical  polarization  can  be 
simplified  as 

{>xE  =  7/oAy  [(i>  •  H)0  -H]^ 

Taking  a  dot  product  on  both  sides  with  z 
z  ■  {i>  X  E)  =  -rjoA"^ H, 

Since  z  X  u  =  —s,  we  have 
-s-E=-r]oA^H, 


Substituting  from  equation  (4)  for  E,  we  get 

.  zxVH 


JU>C 


=  -rioA^H. 
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juerjoA^H,  =  jkoCrA^H, 
i.e.,  the  above  can  be  put  conveniently  in  the  form 


dH, 

du 


-iA:oe.Ar//,  =  0 


IBC  Vertical  Pol. 


A  similar  analysis  for  horizontal  polarization  yields 

IBC  Horizontal  Pol. 


This  may  also  be  obtained  by  resorting  to  duality. 

For  a  perfectly  conducting  material,  we  have  =  0  and 


OH, 

du 

E, 


0  Vertical  Pol. 

0  Horizontal  Pol. 


(16) 


(17) 
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4.  Standard  PE  Derivation 


We  will  make  some  approximations  and  cast  (11)  in  the  form  of  a  parabolic  equa¬ 
tion  which  permits  a  rapid  numerical  solution. 


Let  rp{x,y)  =  — -^uix^y) 


\/x 


Then 


Xl^y  =  Uy  - ,  Uyl  =  U^y  ~  (U^j,  ”  jkoUy)' 


o-jkox 


X  OO 


\/x 


y/x 


^yy  ^y 


7/ 

LL-TT 


(Uji  -  2jkoUr  -  klu)- 


o-jkox 


y/x 


Substituting  into  (17)  we  get 

+  Uyy  -  2jkoUj.  -b  QiUj:  -  2jkoaiU  +  a2Uy  +  (n^  -  l)A-oti  =  0 

or 

Uix  +  fJ'yy  +  (fli  -  2j^o)wi  +  a2Uy  +  -  1  -  k^u  =  0 

If  now  we  impose  the  approximation  that 

|uxi|  ^  -b  4A;q^  |^x|  ) 


we  obtain 


-1 


Ut  = 


(oi  -  2jko) 


Uyy  -b  a2Uy  -b  ^71^  “  1  “  2  j  /Jq  U 


or 


(18) 


This  is  the  exact  form  of  narrow  angle  PE  approximation.  We  would  also  like  to 
express  the  impedance  boundary  condition  in  terms  of  the  ‘u’  functions. 


A 

s 


X  =  s  cos  6  —  i/sin  9 
dx  ..  n  »  - 

x^,  =  —  =  U  ■  Vx  —  U  ■  X 

av 

=  —  sin  0 

k  =  —X  sin  ^  -b  y  cos  0 
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For  vertical  polarization  we  have 


^  -  jkoCrA^H,  =  0  -  jkon^A^H,  =  0 


g-jfcoi 


Hz{x,y)  = 

u{x,y) 

V  X 

dH, 

—  jkox^u 

^-]kox 

dv 

2x3/2 / 

v/? 

^-jkox 

rs^ 

{u^  —  jkox^u) 

X 

■  ■  —  jko{n^A'^  +  x^)u  =  0 

Similarly  for  horizontal  polarization,  we  have 


IBC  Vert.  Pol. 


Uu  -  jko 


IBC  Horz.  Pol. 


We  combine  the  two  by  defining 


'  —  ~  sin  0)  Vert. 

-j f  ^  -  sin  ^  j  Horz. 


and  writing  as 


+  Ciu  =  0 


(19) 


The  parabolic  equation  given  in  (18)  is  valid  for  propagating  angles  close  to  horizon¬ 
tal  (±10°  in  practice)  [1].  To  accomodate  waves  at  higher  angles  we  would  need  a 
wide  angle  parabolic  equation  whose  derivation  is  accomplished  through  a  pseudo¬ 
differential  operator  formalism  [1]. 
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5.  Wide  Angle  PE  Derivation 


Let  US  assume  that  the  media  constants  are  independent  of  horizontal  range  so  that 
a(x,t/)  =  a(?/),  n{x^y)  =  Ti{y).  We  then  have  from  (17) 


Let 


+  u 


yy 


-j2koUi  +  -  ^Uy  +  (n^(j/)  -  l)klu  =  0 
a  ay 


P 


d  \  \  da  d 

dx  ^  \  k^  dy"^  dy  dy 


(Pseudo  differential 
operator  in  y 


Using  this  notation,  the  above  PDE  may  be  expressed  as 

p2-2iA:oP  +  (g"-l)A-2]u=0 

which  we  may  factorize  as 


{P  -  jka  -  jkoQ){P  -jko  +  jkQQ)u  =  0 

To  see  this  we  expand  the  operator  on  the  left  hand  side  to  get 

-  jkoP  -  jkoQP  -  jkoP  -  kl  -  k^Q 

+  jkoPQ  +  Ao'Q'  +  klQ 


(20) 


Since  PQ  =  QP,  we  have  the  desired  result.  The  first  operator  in  (20)  denotes 
an  incoming  wave  (w.r.t.  x)  and  the  second  an  outgoing  wave.  We  retain  only  the 
outgoing  wave  to  obtain 

Pu  = -jkoiQ  -  l)u  (21) 

The  square  root  operator  Q  is  global  in  nature  and  we  would  like  to  make  some 
approximations  to  derive  a  local  operator  from  it.  (It  is  global  because  when  expanded 
in  terms  of  series,  it  will  contain  terms  of  all  derivatives).  Now 


Q  = 


l  l  da  d  I 

kl  a  dy  dy'^  kl  dy^ 


1/2 


Let  us  rewrite  Q  as 


Q  = 


1  +  [^^(j/)  -  i]  +  - 


1  da 


+ 


<C1  normally 


a  d{koy)d{koy)  d{koyy 


small 


1/2 
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which  may  be  expressed  as 

(5  =  vTTv 

where 


V 


—  1  + 


I  da  d 
kla  dy  dy 


ko  dy'^ 


a  small  operator! 


Treating  V  as  an  algebraic  factor  <  1,  we  may  derive  the  following  rational  approxi¬ 
mation  {pade{\,  1)) 


Q=  vTTV 


SO  that 


i+ii/ 

4  +  3T 
4  +  T 

2V 


(Claerbout) 


=  1  + 


2V 
4  +  V 


Substituting  into  (21)  we  arrive  at 

Pu  =  jko{Q  -  l)u 


—2jkoV 

Pu  =  — - —u 

4  +  V 


or 


Using  the  fact  that  a2  —  — 

a 


(4  +  V)Pu  =  —2jkQVu 
da 

— ,  we  get 
dy 


/  2  o\  I  2  ^ 

(n  + 


Ui  =  -2jko 


{n^  —  1)^0  d"  ^2"^ — I" 


dy  dy'^ 


(22) 


The  above  equation  is  a  wide-angle  parabolic  equation  valid  for  propagation  up  to 
±20°  [1].  Other  approximations  could  be  obtained  by  considering  higher  order  pade 
approximants. 
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6.  Boundary  Condition  on  the  Upper  Boundary 


To  truncate  the  computational  domain,  we  consider  a  point  high  enough  where  the 
atmosphere  is  homogeneous  with  n  =  1.  The  governing  equation  in  the  homogeneous 
region  becomes 


2ko 

Let  us  derive  boundary  conditions  on  a  horizontal  interface  y  =  yo-  For  the  sake 
of  simplicity  we  will  derive  boundary  conditions  of  the  mixed  type  on  the  interface 
y  =  0  (instead  of  y  =  yo)  given  the  initial  data  on  the  line  x  =  0  and  boundary  data 
on  the  line  y  =  0.  Our  derivation  is  based  on  the  use  of  Fourier  sine  transforms  as 
suggested  in  [2].  Although  the  basic  philosophy  of  our  approach  coincides  with  that 
in  [4],  some  of  the  details  and  the  final  results  are  slightly  different  from  the  latter. 


Consider  the  parabolic  equation  Uyy  —  2jkuj:  =  0  in  a  homogeneous  region  x  >  0, 
y  >  0,  where  A:  is  a  complex  constant, 


subject  to  the  initial  condition  u(0,  y)  =  /(y),  0  <  y  <  cx),  and  the  boundary  condition 
u(x,0)  =  y(x),  0  <  X  <  oo.  We  assume  that  u(x,oo)  — >  0  and  Uy(x,oo)  — >  0.  The 
equation 

Uyy  =  2jkUj:  (23) 

is  of  Schrodinger’s  type.  We  will  treat  the  lossless  case  having  a  real  value  of  k  as 
the  limiting  case  of  the  lossy  problem  having  k  =  ko  —  jt^  e  >  0.  We  will  solve  the 
problem  using  Fourier  sine  integral. 


Let 


[2 

Us{x,\)  =  \-  u{x,y)s\nXydy 
V  TT  Jo 


(24) 
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Then  using  integration  by  parts,  we  see  that 


TT  Jo 


Uyy{x,  y)  sm  Xy  dy  =  \  -  luy{x,y)  sin  Xy 


TT 


y=0 


—Xu{x,  y)  cos  Xy 


y=0 


r 

—  y  u{x ^  y)  sin  Xy  dy 


Because  Uy{x,  cxd)  — >  0  and  u{x,  oo)  — >  0,  we  have 
(2  f°o 


—  Uyy{x,y)  sin  Xu  dy  =  1 /— Au(a:,  0)  —  A^[/j(x,  A) 


TT  Jo 


TT 


=  \I-Xg{x)  -  X'^Us{x,X) 


(25) 


Multiplying  both  sides  of  (23)  with  y2/7rsinAy,  integrating  over  3/  =  0  to  00  and 
making  use  of  (25),  we  have 


or 


^^Xg{x)  -  X'^Usix,  A)  =  2jk-^Us{x,  A) 


Q^Us{x,X)  +  ^jf^U,{x,X)  - 


We  may  rewrite  the  above  equation  as 


d 


A  2 


-gix)e 


{X^/2jk)T 


(26) 


2jk  V  tt" 

Replacing  the  dummy  variable  x  in  (26)  with  r  and  integrating  both  sides  over  r  =  0 
to  X,  we  arrive  at 


A  ^ 


^^0  2;/:  V  TT  Jr=o 


g(T)e^^y^jk)rdT 


or 


But 


Us{x,X)  =  Us{0,X)e-^^  + 


-(\'^l2]k)x  ^  /2 


2jk 


V  TT  Jt=Q 


UsiO,X)  =  \  -  f  u{0,y)  sin  Xydy 

V  TT  Jo 

=  r  f(y)im\ydy  =  F.(X) 
\  TT  Jo 
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U,{x,\)  =  ^  A  /2  r  (27) 

ZjK  \  TT  Jr=0 

Finally,  taking  the  inverse  sine  transform  on  (27)  we  get 

u{x,y)  =  W-  f  sm  XyUs{x,X)dX 

V  TT  Jo 

V  TT  Jx=0 

+  -7rT  r  A  sin  Ay  T  dA 

7r2jA:JA=o  Jt=o 


=  —  /  /(r)  sin  At  sin  Ay  (7A 

TT  7a=0  7t=0 

+  -^  r  y(r)  /°°  AsinAye<^'/2^'=)<^-^)<fAdT 

TTZjkJrzzO  JX=0 

1  /‘OO  /*oo 

=  —  /(r)  /  [cos(r  —  2/)A  —  cos(r  +  j/)A]  (fr 

TT  Jt=o  Ja=o 

+  -  r  5(^)4  f-|-l  /°°  cosAye(^'/2^*^)(^-")c7A^fr 

TT  Jt=0  V  ^2//  ^^^=0 

=  —  /(■’■)/  [cos(y  —  t)A  —  cos(y  +  t)A]  dr 

X  7t=0  ja=o 


Defining 


/\  (x,y;  xo,yo)  =  /  cos(y  -  yo)Ae  (^V2jfc)(x  for  x  >  lo, 

^0 

we  write  the  expression  for  u(x,y)  as 

1  r®® 

u(x,y)  =  -  /(r)  (/\(x,y;  0,t)  - /\(x,y;  0, -T)]dT 

TT  Jt=0 

1  5 

TryK  Jt=o  ay 

We  now  evaluate  the  integral  for  K.  Consider 

I{a)  =  f  ,  x  >  xq 

Jx=o 


(28) 
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Jx=o 

Jx^o 


Because  of  the  exponential  decay,  we  may  differentiate  under  the  integral  sign  to 
obtain 

A/(o)  =  -  /°° 
da  Jx=o 

g-(A2(i-ro)/2|fc|^)(e-ifco)l 


=  /  si 

Jx=o 


sin  aA 


_a_ 

aA 


{x  -  io)(e  -jko)/\k\‘^ 


d\ 


sin 


a|^P 


(x  -  Xo){e  -  jko)/\k\'^ 

r  cos 
Jo 

ja\k\‘^ 


=0  {x  -  Xo){e-  jko) 


{x  —  Xo)k 


-■L 


{x  -  Xo)' 


dlia)  jak  ,,  _ 


or 


i 


,jo^<:/(2(i-io)) 


=  0  =>  I{a)  =  /(a  =  0)e--’“'''/(2(x-xo)) 


Now  /(a  =  0)  can  be  written  as 

/(a  =  0)  =  r 
io 

Let  US  view  this  integral  in  the  complex  A-plane. 


(29) 


I 


For  the  integral  to  converge  for  {x  —  xq)  >  0,  Re[X^{e.  —  j/:o)]  >  0,  i.e., 

Re[X^-jk-)]  >0=>  Im{X^k’)  >  0 


or 


0  <  Arg  {X^k*)  <  TT 


or 


or 


Now 


0  <  2  Arg  (A)  —  Arg  (k)  <  ir 


^  Arg  (k)  <  Arg  (A)  <  ^  ^  Arg  {k) 


Arg  [k)  =  —  tan 


-1 


^0 


ko  ko 


<  1 


/(a  =  0)=  / 

Jc 

where  C  lies  in  the  region  of  convergence. 

In  particular,  let  us  choose  a  line  from  0  to  oc  along  the  line  Cuj  defined  by 


(30) 


On  this  path, 


.  'X  1  ...  TT  6 

Ar6(A)  =  j  +  -Arg(i-)  =  j-5j^ 


Ttjk 


2(x  -  Xo) 
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/(a)  = 

V  2{x  —  Xo) 


(32) 


Substituting  this  into  (28),  the  field  u{x,y)  for  the  initial-boundary  value  problem  is 
given  by 

u{x,y)  =  i  /“ 

TT  Jt=0  V  2x  ^  ^ 


._L  r  glr)±iJ 

Trjk  Jtso  dy  2{i  —  t)  j 


It  is  easy  to  see  that 

d^K 


dy^ 


=  r  -A^  cos[(y  - 
Jo 


and,  so  for  x  ^  xq 


-n  1  d^K 

dy^  ^  dx  dx  2jk  dy^ 


which  is  the  same  equation  satisfied  by  u.  Now 

/°°  /(r)  _  g-iA(v+r)V(2r) 

V  27rx  Jt=o 

1  d 

- rr  /  9{T)-^I<{x,y,  T,0)dT 

7r?fcJr=o  uy 


Ijk  r 

y— 0+ 

V  27rx  do 

r  l(r)  — 

Jo  X 


.g-jArV(2x)^^ 


Jo  ar 


TTJ 


cfr 


y-^o+ 


27rx 


irjk 


71 


dr 


y-»0+ 


(33) 


(34) 
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in  view  of  (34).  Furthermore,  we  note  that 


2/5  yo)  .^'^(^02/5^02/0) 

Using  this  and  evaluating  the  first  integral  by  parts,  we  get 
d 


dy 


u(x,y) 


y— 0+ 


+-  /  y(T)^A'(x,y;  r,0)|  dr 

71"  Jt=0  OT  . 


y^0+ 


2 

+- 

7r 


i?(r)A:(x,y;T,0) 

,_o  -  /  yT(7-)/\(x,0;  r,0)dr 

V-0+  ''^=° 

=  +J—m  + 

V  TTX 


I2jk 


TTX 


Jo 


-^S(O), 


hjk  _  grjr) 

V  TT  Jt=0  \/x  —  r 


2x 


v/a 


dr 


I2jk 


TTX 


1/(0) -5(0)1 


^  /t(t),-,i„=/|,. 

TT  1  Jo 


Jo  vx  —  r 


dr 


\/x  Jo  y/a 

From  the  compatibility  conditions  on  the  initial  and  boundary  values,  we  have 


Therefore 

du 


dy 


{x,y) 


limu(x,0)  =  5r(0)  =  limu(0,y)  =  /(O) 

X— *-0  y— *^0 


^  r  r  -plLdr 

y/X  Jt=0  Jt=0  y/x  —  T 


y-^0+ 


9r{T) 


It  is  easy  to  note  from  (35)  that  for  k  =  ko  —  jt, 


du 


=  0. 


y^0+ 


(35) 
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No  matter  how  small  the  c  is,  we  will  use  the  same  equality  for  the  limiting  case 
of  e  — >  0.  We  will  evaluate  the  integrals  above  approximately  by  replacing  the 
derivatives  with  differences. 


y  t 


Ay 


I 


u. 


u. 


u. 


u 


\ 


X 


Consider  initial  data  on  the  line  i  =  0.  Let  us  assume  that  this  initial  data  is  known 
on  a  uniform  grid  Pm  =  mAy,  m  =  0, 1,  •  •  •.  We  approximate  the  derivative  in  the 
interval  (ym-iiym)  by  the  forward  difference  formula 


dy  Ay 


y  e  (2/m-l,2/m) 


where  =  w(0,y,n)- 


Then 


Let 
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We  then  have 


yJkl(Tx)y 

m 


= 


yyk/{TTx)ym-l 

y/k/(iTx)y  m 
yjkl(xx)y  m  — 1 


where  F(/i)  =  complex  Fresnel  integral  [8] 


1 

-r 

yJX  JrzzO 


E( ^m— l\ 

A  Ay  jV/^ 


J/n 


2/m  —  1 


(36) 


Consider  now  the  boundary  data  on  the  line  j/  =  0.  Assume  that  we  have  a  non- 
uniform  grid  0, xi ,  j:2,  *  •  • , xa/  =  x.  On  the  interval  we  approximate  the 

derivative  as 

9r{T)  = 

^m— 1 

where  u’”  =  u(xm,0).  At  the  origin  we  have  u°  =  Uq. 


_  ym-l 


We  evaluate  the  second  integral  in  (35)  as 

5,(r)  1 


r  J^<lr  =  £  " .  r  -7^<lr 

Jo  y/ X  T  m=l  ^m  — 1  J^m—\  y/^  ^ 


=  i: 


Xm  ^m-1 


(^—2\/x  ~  ) 


3^171  —•  1 


=  2i: 


u”*  -  u"*"* 


—  1  —  1 


m  =  l 


M 


^a/x  im  — 1  y/^  ) 


=  2E 


u'"  -  u'"“* 


—  2  \/^ Af  ^m  — 1  d”  \/^A/ 


(37) 


Using  (36)  and  (37)  in  (35)  we  get 


y=0+ 


f^Jk  [tF  Um  ^m-l  i  Tp  (  I  ^ 


-?/m  -  F 


TTX  Af 


‘Vm  —  1 


M 


-2E 


u'"  - 


—2  A/  ^771  d”  \/^Af  —  1 
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Extracting  out  the  m  =  M  term  from  the  latter  and  writing 


we  have 


■^0+ 


du^ 

dy 


du^  ^  8jk 
dy  V  TT  y/XM  -  XM-1 


1 


F 


-F 


_  ym-\ 


_ 

^  m=l  \/^M  “t*  \/^M~~~^rn— T 


(38) 


The  above  equation  is  the  discrete  version  of  a  continuous  boundary  condition  of  the 


form 

du  ,  ,  ,  , 

(39) 

—  -t-  r(x)u  =  s(x) 
dy 

where 

and 

r(x  =  \  - 

]l  TT  y/x  -  XM-l 

(40) 

m  =  l  ^ 


F\yl-v^\-F 


k 

J/m  — 1 


M-1 


-\/^E 


u*”  —  u”*  ^ 


TT  rn—1  *^m  “i”  ^m  — 1 


and 


F(X)  = 

Jo 

is  the  complex  Fresnel  integral  [8]. 


(41) 


(42) 


For  efficient  implementation,  the  PDE  and  the  various  boundary  conditions  will  be 
transformed  into  a  curvilinear  coordinate  system  generated  by  setting  the  lower  ir¬ 
regular  boundary  as  a  constant  curve  curvilinear  coordinate. 


23 


( 


7.  Transformations  to  a  Curvilinear  Coordinate  System 


Consider  the  narrow  angle  parabolic  equation  with  range  dependent  refractive  index 
(oi  ^  0)  given  in  (18) 


1  [  .  d 

together  with  the  tropospheric  boundary  condition 

0 

Qy  »  J/o)  d”  )^(^m  )  J/o)  •®(^m  )  2/o) 

and  the  impedance  boundary  condition  on  the  irregular  boundary 


*7 - 7^ 


+  ciu  =  0 


0^/^OL 


^0\A/VV,cC^Ay 


We  will  transform  these  to  a  curvilinear  coordinate  system  {(,ri) 
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We  assume  that  we  have  a  transformation  of  the  following  form 

The  various  metrics  needed  in  the  transformed  equations  are  [5] 

9n  =  +  ,  gi2  =  +  y^Vr, . 

We  then  have  with  y/g  =  ^  0 

Ux  =  -^iyr,U(  -  y^Ur,)  = 

y/9  ^(Vr) 

Uy  =  \i{-XnU(_-\-  X(_Ur))  =  — 

y/9  Vr, 


1  d 


^yy  ~ 


(“"'I 

1 

0.2 

^r}T}  ^77 

L  Vr)  J 

Substituting  into  the  PDE  we  get 

{u^Vv  -  ^vva  =  tt- 


^(.Vn 


{2jko  -  Ci) 


«  I  Vriv  I  ^vv  \ 

ijix  +  a2 - :rUr,  d - j-  / 

Vr,  yl  J 


or 


XirOi 

- -U  + 


{2jko  —  ai) 


®2  yrjTt 


+ 


.yv  J/^  /  (2i^o -fli)  yr, 


Xf 

"h  fr,  •  i,  ^  2^'^^ 

{2jko  -  ai)y2 


Letting 


6i  = 


{2jko  —  a-i) 


6,  = 


yv 


63  — 


02  - 


X^ 


y^v 


1 


y^  J  (2iA:o  -  Oi) 


+  y? 


{2jko  -  ai)y2 


we  express  the  PDE  as 


—  6]  u  -f-  62^77  “f“  ^3^ 


777? 


(43) 


The  normal  derivative  on  a  y  =  constant  line  is 


r  I  \  3u  912 

u^(T]  =  const.)  =  ./ — u,j - ix£ 

V  9  y/99n 


\J^\  +  y\ 


^iyv-yi^n  ”  ^Jx\  +  y|(x^y^  -y^x^) 


+  y^y») 


(44) 
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Variation  of  an  arbitrary  vector  r  along  the  t)  =  const,  coordinate  is 

=  x^x  +  y^y  = 

The  unit  vector  along  the  tangent  on  r;  =  const,  is  then 


s  = 


x^x  y^y 


Vi 


Defining 


cos  6  = 


sin  6  = 


x^ 


\AF+ 


Vi 


Vi 


we  express  (44)  as 


til/  = 


{xr,  COS  0  +  yr,  sin  6)u(^ 


yr,  cos  6  -  Xr,s\n  9  +  y|(y^  COS  0  —  Xr^  sin  0) 


With  these  substitutions,  the  boundary  condition  at  the  bottom  boundary  h 
0  gets  transformed  to 

(x^cos0  +  y^sin0) 

- , . .  . . +  [yr)  cos  U  —  Xrj  Sin  U)C\  U  =  (J  on  7]  =  U 


\lA  + 


Vi 


For  =  0  and  using  y^x^  +  yj  =  x^/  cos  6,  we  get 


V 

u-n - ^  sin  0  cos  4- Cjj/^  cos  =  0  @  rj  =  0 


Finally  at  the  top  boundary,  we  have 

^v/Vr)  +  ru  =  s 


or 


Ur)(^m)  A  )  T  r(^,„,  A^)yr;(^m)^(^m5  -V) 


1/  +  C]  U  — 

(45) 


(46) 


(47) 
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8.  Generation  of  Curvilinear  Coordinate  System 


Consider  a  piece-wise  linear  ground  profile  and  a  horizontal  upper  boundary 


rjtiN 


■  -  ...  ^ 

^  ^  £ 

f— ^ 1 

B' 

^  t 

c' 

D 

^ 

-J - 

Arj 

tl 

7 - —r* 

^ 

-7 - 7—^^ 

I  7  ■  ^ 

> - y - T  ^ - 7 - /  -  •  /  r  "  ^ 

A  6  C  D  ' 


Non-Rectangular, 
Uniform  Mesh 


Rectangular, 
Uniform  Mesh 


Physical  Domain 


Computational  Domain 


^ _ 6 _ 6 — £■ — « 


N 


We  generate  the  {x,y)  coordinates  of  an  interior  point  by  using  a  linear  interpolation. 
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Letting  Ax,  =  (x.+i  - 1,),  Ay,  =  (y,+i  -  j/,),  and  A^,  =  (^.+1  -  we  write 


x  =  x,  +  -  (,) 


y 


<^ywc  c\ 


+  jv^« 


(.<(<  f.+l 

0  <y  <N 


At  any  interior  point,  (i  <  (  <  6+11  the  various  metrics  are  evaluated  as 
Ax, 


x^  = 


=  1 


Ae.  ’ 


X,,  =  0 


nJ 


n  \  Aj/, 


A6  ’  N 


Ay,  f  ^  ^  ^ 

yo  -  y.  -  -  ix) 


(48) 


(49) 


9  V'n'n  —  0 


At  the  boundary  points,  we  use  the  central  difference  formulas  to  arrive  at 


x,+2  -  a:,  Ax,+1  +  Ax,- 

(50) 

6+2  —  6  Aif,+i  +  A^, 

II 

II 

A  Ay.+i  +  Ay, 

V  n)a(,+,  +  a(. 

(51) 

Note  that  the  analytical  expressions  yield  a  discontinuous  value  for  these  derivatives 
at  the  boundary  points.  We  use  the  analytical  expressions  only  to  generate  grid 
points  and  use  the  central  difference  formulas  to  arrive  at  the  derivatives  w.r.t. 

In  other  words,  once  the  grid  points  are  generated  on  the  lines  AA\  BB\  . . . ,  etc., 
we  assume  that  the  space  is  smoothly  connected  through  the  grid  points.  In  the 
numerical  implementation  using  Crank-Nicolson  implicit  scheme  [6],  the  metrics  are 


needed  at  the  midpoint  w.r.t.  i.e.,  at  ^  =  6  +  A672  and  the  interior  point  formulas 
are  applicable.  For  a  uniform  mesh  in  the  computational  domain,  A^,  =  1,  rj  =  q, 
q  =  0,1,2,---,A^ 

x^  = 

Ax, 

(52) 

y(.  = 

(i  yv)Ay.  =>  y^(y  +  l)  ydy)= 

(53) 

yr,  = 

^  f..  yi+i  +  yA 

NV^  2  ) 

(54) 

y  = 

(,  fyx+i+yx\  ,  q 

['  N>  \  2  j  +  A'*'” 

= 

2  -^m  =  yo  +  (y  ^)s'i 

(55) 

so  that  y{q  +  1)  -  yiq) 

=  Vv' 
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9.  Numerical  Implementation  by  Crank-Nicolson  Scheme 

Consider  the  narrow  angle  PE  together  with  the  boundary  condition 

=  biU  +  b2Ur,  +  bzUr^r,  (PDE) 

—  I  —  sin  0  cos  0  )  Uf  +  c\yrf  cos  Ou  =  0  at  77  =  0  (BC2) 

/ 

Ur,  +  ry^u  =  Pr,  s  at  7;  =  iV  (BCi) 

with 

Ur,(0,N)=0 

We  would  like  to  implement  the  above  using  a  Crank-Nicolson  implicit  scheme  [6] 


i 

i-' 


We  use  the  notation  u^  =  u{(p,Tjg)  =  u{xp,yg) 

The  various  derivatives  assuming  =  Arj  =  1  are 

(^P-l/2)  ~ 

=  7  [<+l  -  <-l  +  <+i  -  <-l] 


or 


p— 1  I 

^9+1  + 


u 


9+1 


p-i 

g-l 


+  <-i 


)1 
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Uvri  K+1  “  -  2uP  ^  +  UP_}] 

Substituting  into  the  PDE,  we  get 

1  /  1  \  u?  + 

<-<  =  '•'V’-r^)— — 


(<+l  +  <+i  -  <-l  - 


1 


«  ^  (^9+1  2u^  +  ti^_i  +  2u^ 


Rearranging  the  terms  we  get 

1  (b 


bx 


1 


\  ( bo 


bi' 


+^l^  +  ^3|<;l  +  |i-63  +  ^|<-’  +  iU3-4  <:}  =  o 


2  V  2 


Now  let 


a 


1  N^P-1/2 

^3+2^2 


0  =  (^3-1 


9 

P-1/2 


7  = 


P-1/2 


We  then  have  for  p  =  1 , 2,  •  •  ■,  q  =  1, 2.  •  •  •,  A'^  —  1 


0!U„ 


-  (1  +  ^)<  +  7<-i  =  -  (1  -  -  7<:i 


,p-i 


,p-i 


(56) 


However,  we  will  extend  the  applicability  of  this  equation  over  q  =  0, 1,--  A^  to 
accomodate  the  derivative  boundary  conditions  on  the  lower  and  upper  boundaries. 
For  g  =  0,  we  have 

aul  -  (1  +  )3)ug  +  7u^i  =  -au?“'  -  (1  -  ^)uo~^  -  7U-V  (57) 

for  g  =  A^,  we  have 

-  (1  +  /?)u^  +  7U%_■^  =  -  (1  -  -  ju’l^-x  (58) 
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From  the  BC2  at  the  lower  boundary  we  have 


-! - f- 


ft--' 


7 - 7^ — X — P - y - 7-  ^  ~  0 

‘  1 

-1 -  1)---' 


I 

I 


p-1  ^,P  ■ 


.  .  -  .  r  ^  /  \  P’’l/2 

u\  -i/li  ,  Vv  ^  n  n\  (  V  p-i\ 

- ^ - 1- - ^ sm^cos^  I  ) 


+(c,„cos«r'^f^^^^ 


= 0 


Multiplying  this  with  4(7)0  adding  to  (57) 

(a  +  'y)u\  -  ( 1  +  /?  -  27 


y 

Ciy„  cos  6  —  2—  sin  0  cos  9 


Un 


=  -(q  +  7)ui  '  -  I  1  -  ^  +  27 


Letting 


o'  =  0  +  7 


V 

Ciyjj  cos  6  +  2—  sin  ^ cos  9 


u: 


p-i 


13'  =  13  -  2'yyr,  cos  0  |  Ci  -  j 


/3"  =  I3  -  2'ry^  cos  9  ici  + 


2  sin  0' 


we  may  write  the  above  equation  as 

a'uj  -  (1  +  ^')ul  =  -a'up*  -  (1  -  I3")ul~^ 

From  the  BCi  at  the  top  boundary  we  have 

- i - * -  N+l 


X  L. 


X  '  ^ ^  N 


|>-l 


N-l 


(59) 
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u 


P-1 

A^+1 


—  U 


P-1 

A^-l 


+ 


^A^-H  ~  ^ 


A^-l 


+ 


pP  -f  pP  *  u%  4- 


■Vv' 


N  SP  +  S 

-  =  Vr,— 


P-1 


Multiplying  this  with  —4(0-)^  and  adding  to  (58),  we  get 

-  (l  +  /S  +  [pP  +  rP"*]y,,a)  u%  +  +  4a)u%_^ 

=  -  (l  -  /0  -  [pP  +  pP"*]j/,,q)  -  (7  +  4o)u^"_\  -  2y„(s'’  +  sP~*) 

Letting 

A  =  ,9  +  r/,,a[p'’ +  pP“’] 

7^  =  7  +  4a 

we  write  the  above  equation  as 

(1  +  A)uJ,  -  7'uJ,.,  =  (1  -  A)uJ,“'  +  7'u5r-'i  +  2!/„('''  +  »"’■') 


(60) 


We  augment  the  equations  given  in  (56)  for  9  =  1 , 2, . . . ,  A'^  —  1  with  (59)  and  (60)  for 
q  =  0  and  N,  respectively  to  define  for  g  =  0, 1,2, . . . ,  A^.  The  system  of  equations  so 
defined  can  be  expressed  as  a  matrix  equation  of  the  form 


A  A  A 
A  A  A 


A  A 


- 

ul 

r 

,  p 

Ul 

,  p 

- 

- 

A  A 
A  A  A 
A  A  A 


A  A 


r  ,  p-i  1 

p—  1 

,  P-1 

(61) 


where  X  denotes  a  non-zero  entry.  The  tridiagonal  matrix  on  the  left  hand  side  of 
(61)  can  be  inverted  efficiently  to  yield  a  solution  on  line  ^  in  terms  of  the  field 
values  on  the  line  (  =  (fp-i-  Equation  (61)  can  be  used  to  march  forward  in  range 
starting  from  initial  data  specified  on  =  0. 
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